14a238c70SJohn Marino /* mpfr_sin_cos -- sine and cosine of a floating-point number
24a238c70SJohn Marino
3*ab6d115fSJohn Marino Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4*ab6d115fSJohn Marino Contributed by the AriC and Caramel projects, INRIA.
54a238c70SJohn Marino
64a238c70SJohn Marino This file is part of the GNU MPFR Library.
74a238c70SJohn Marino
84a238c70SJohn Marino The GNU MPFR Library is free software; you can redistribute it and/or modify
94a238c70SJohn Marino it under the terms of the GNU Lesser General Public License as published by
104a238c70SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
114a238c70SJohn Marino option) any later version.
124a238c70SJohn Marino
134a238c70SJohn Marino The GNU MPFR Library is distributed in the hope that it will be useful, but
144a238c70SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
154a238c70SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
164a238c70SJohn Marino License for more details.
174a238c70SJohn Marino
184a238c70SJohn Marino You should have received a copy of the GNU Lesser General Public License
194a238c70SJohn Marino along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
204a238c70SJohn Marino http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
214a238c70SJohn Marino 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
224a238c70SJohn Marino
234a238c70SJohn Marino #define MPFR_NEED_LONGLONG_H
244a238c70SJohn Marino #include "mpfr-impl.h"
254a238c70SJohn Marino
264a238c70SJohn Marino /* (y, z) <- (sin(x), cos(x)), return value is 0 iff both results are exact
274a238c70SJohn Marino ie, iff x = 0 */
284a238c70SJohn Marino int
mpfr_sin_cos(mpfr_ptr y,mpfr_ptr z,mpfr_srcptr x,mpfr_rnd_t rnd_mode)294a238c70SJohn Marino mpfr_sin_cos (mpfr_ptr y, mpfr_ptr z, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
304a238c70SJohn Marino {
314a238c70SJohn Marino mpfr_prec_t prec, m;
324a238c70SJohn Marino int neg, reduce;
334a238c70SJohn Marino mpfr_t c, xr;
344a238c70SJohn Marino mpfr_srcptr xx;
354a238c70SJohn Marino mpfr_exp_t err, expx;
364a238c70SJohn Marino int inexy, inexz;
374a238c70SJohn Marino MPFR_ZIV_DECL (loop);
384a238c70SJohn Marino MPFR_SAVE_EXPO_DECL (expo);
394a238c70SJohn Marino
404a238c70SJohn Marino MPFR_ASSERTN (y != z);
414a238c70SJohn Marino
424a238c70SJohn Marino if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
434a238c70SJohn Marino {
444a238c70SJohn Marino if (MPFR_IS_NAN(x) || MPFR_IS_INF(x))
454a238c70SJohn Marino {
464a238c70SJohn Marino MPFR_SET_NAN (y);
474a238c70SJohn Marino MPFR_SET_NAN (z);
484a238c70SJohn Marino MPFR_RET_NAN;
494a238c70SJohn Marino }
504a238c70SJohn Marino else /* x is zero */
514a238c70SJohn Marino {
524a238c70SJohn Marino MPFR_ASSERTD (MPFR_IS_ZERO (x));
534a238c70SJohn Marino MPFR_SET_ZERO (y);
544a238c70SJohn Marino MPFR_SET_SAME_SIGN (y, x);
554a238c70SJohn Marino /* y = 0, thus exact, but z is inexact in case of underflow
564a238c70SJohn Marino or overflow */
574a238c70SJohn Marino inexy = 0; /* y is exact */
584a238c70SJohn Marino inexz = mpfr_set_ui (z, 1, rnd_mode);
594a238c70SJohn Marino return INEX(inexy,inexz);
604a238c70SJohn Marino }
614a238c70SJohn Marino }
624a238c70SJohn Marino
634a238c70SJohn Marino MPFR_LOG_FUNC
644a238c70SJohn Marino (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
654a238c70SJohn Marino ("sin[%Pu]=%.*Rg cos[%Pu]=%.*Rg", mpfr_get_prec(y), mpfr_log_prec, y,
664a238c70SJohn Marino mpfr_get_prec (z), mpfr_log_prec, z));
674a238c70SJohn Marino
684a238c70SJohn Marino MPFR_SAVE_EXPO_MARK (expo);
694a238c70SJohn Marino
704a238c70SJohn Marino prec = MAX (MPFR_PREC (y), MPFR_PREC (z));
714a238c70SJohn Marino m = prec + MPFR_INT_CEIL_LOG2 (prec) + 13;
724a238c70SJohn Marino expx = MPFR_GET_EXP (x);
734a238c70SJohn Marino
744a238c70SJohn Marino /* When x is close to 0, say 2^(-k), then there is a cancellation of about
754a238c70SJohn Marino 2k bits in 1-cos(x)^2. FIXME: in that case, it would be more efficient
764a238c70SJohn Marino to compute sin(x) directly. VL: This is partly done by using
774a238c70SJohn Marino MPFR_FAST_COMPUTE_IF_SMALL_INPUT from the mpfr_sin and mpfr_cos
784a238c70SJohn Marino functions. Moreover, any overflow on m is avoided. */
794a238c70SJohn Marino if (expx < 0)
804a238c70SJohn Marino {
814a238c70SJohn Marino /* Warning: in case y = x, and the first call to
824a238c70SJohn Marino MPFR_FAST_COMPUTE_IF_SMALL_INPUT succeeds but the second fails,
834a238c70SJohn Marino we will have clobbered the original value of x.
844a238c70SJohn Marino The workaround is to first compute z = cos(x) in that case, since
854a238c70SJohn Marino y and z are different. */
864a238c70SJohn Marino if (y != x)
874a238c70SJohn Marino /* y and x differ, thus we can safely try to compute y first */
884a238c70SJohn Marino {
894a238c70SJohn Marino MPFR_FAST_COMPUTE_IF_SMALL_INPUT (
904a238c70SJohn Marino y, x, -2 * expx, 2, 0, rnd_mode,
914a238c70SJohn Marino { inexy = _inexact;
924a238c70SJohn Marino goto small_input; });
934a238c70SJohn Marino if (0)
944a238c70SJohn Marino {
954a238c70SJohn Marino small_input:
964a238c70SJohn Marino /* we can go here only if we can round sin(x) */
974a238c70SJohn Marino MPFR_FAST_COMPUTE_IF_SMALL_INPUT (
984a238c70SJohn Marino z, __gmpfr_one, -2 * expx, 1, 0, rnd_mode,
994a238c70SJohn Marino { inexz = _inexact;
1004a238c70SJohn Marino MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
1014a238c70SJohn Marino goto end; });
1024a238c70SJohn Marino }
1034a238c70SJohn Marino
1044a238c70SJohn Marino /* if we go here, one of the two MPFR_FAST_COMPUTE_IF_SMALL_INPUT
1054a238c70SJohn Marino calls failed */
1064a238c70SJohn Marino }
1074a238c70SJohn Marino else /* y and x are the same variable: try to compute z first, which
1084a238c70SJohn Marino necessarily differs */
1094a238c70SJohn Marino {
1104a238c70SJohn Marino MPFR_FAST_COMPUTE_IF_SMALL_INPUT (
1114a238c70SJohn Marino z, __gmpfr_one, -2 * expx, 1, 0, rnd_mode,
1124a238c70SJohn Marino { inexz = _inexact;
1134a238c70SJohn Marino goto small_input2; });
1144a238c70SJohn Marino if (0)
1154a238c70SJohn Marino {
1164a238c70SJohn Marino small_input2:
1174a238c70SJohn Marino /* we can go here only if we can round cos(x) */
1184a238c70SJohn Marino MPFR_FAST_COMPUTE_IF_SMALL_INPUT (
1194a238c70SJohn Marino y, x, -2 * expx, 2, 0, rnd_mode,
1204a238c70SJohn Marino { inexy = _inexact;
1214a238c70SJohn Marino MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
1224a238c70SJohn Marino goto end; });
1234a238c70SJohn Marino }
1244a238c70SJohn Marino }
1254a238c70SJohn Marino m += 2 * (-expx);
1264a238c70SJohn Marino }
1274a238c70SJohn Marino
1284a238c70SJohn Marino if (prec >= MPFR_SINCOS_THRESHOLD)
1294a238c70SJohn Marino {
1304a238c70SJohn Marino MPFR_SAVE_EXPO_FREE (expo);
1314a238c70SJohn Marino return mpfr_sincos_fast (y, z, x, rnd_mode);
1324a238c70SJohn Marino }
1334a238c70SJohn Marino
1344a238c70SJohn Marino mpfr_init (c);
1354a238c70SJohn Marino mpfr_init (xr);
1364a238c70SJohn Marino
1374a238c70SJohn Marino MPFR_ZIV_INIT (loop, m);
1384a238c70SJohn Marino for (;;)
1394a238c70SJohn Marino {
1404a238c70SJohn Marino /* the following is copied from sin.c */
1414a238c70SJohn Marino if (expx >= 2) /* reduce the argument */
1424a238c70SJohn Marino {
1434a238c70SJohn Marino reduce = 1;
1444a238c70SJohn Marino mpfr_set_prec (c, expx + m - 1);
1454a238c70SJohn Marino mpfr_set_prec (xr, m);
1464a238c70SJohn Marino mpfr_const_pi (c, MPFR_RNDN);
1474a238c70SJohn Marino mpfr_mul_2ui (c, c, 1, MPFR_RNDN);
1484a238c70SJohn Marino mpfr_remainder (xr, x, c, MPFR_RNDN);
1494a238c70SJohn Marino mpfr_div_2ui (c, c, 1, MPFR_RNDN);
1504a238c70SJohn Marino if (MPFR_SIGN (xr) > 0)
1514a238c70SJohn Marino mpfr_sub (c, c, xr, MPFR_RNDZ);
1524a238c70SJohn Marino else
1534a238c70SJohn Marino mpfr_add (c, c, xr, MPFR_RNDZ);
1544a238c70SJohn Marino if (MPFR_IS_ZERO(xr)
1554a238c70SJohn Marino || MPFR_EXP(xr) < (mpfr_exp_t) 3 - (mpfr_exp_t) m
1564a238c70SJohn Marino || MPFR_EXP(c) < (mpfr_exp_t) 3 - (mpfr_exp_t) m)
1574a238c70SJohn Marino goto next_step;
1584a238c70SJohn Marino xx = xr;
1594a238c70SJohn Marino }
1604a238c70SJohn Marino else /* the input argument is already reduced */
1614a238c70SJohn Marino {
1624a238c70SJohn Marino reduce = 0;
1634a238c70SJohn Marino xx = x;
1644a238c70SJohn Marino }
1654a238c70SJohn Marino
1664a238c70SJohn Marino neg = MPFR_IS_NEG (xx); /* gives sign of sin(x) */
1674a238c70SJohn Marino mpfr_set_prec (c, m);
1684a238c70SJohn Marino mpfr_cos (c, xx, MPFR_RNDZ);
1694a238c70SJohn Marino /* If no argument reduction was performed, the error is at most ulp(c),
1704a238c70SJohn Marino otherwise it is at most ulp(c) + 2^(2-m). Since |c| < 1, we have
1714a238c70SJohn Marino ulp(c) <= 2^(-m), thus the error is bounded by 2^(3-m) in that later
1724a238c70SJohn Marino case. */
1734a238c70SJohn Marino if (reduce == 0)
1744a238c70SJohn Marino err = m;
1754a238c70SJohn Marino else
1764a238c70SJohn Marino err = MPFR_GET_EXP (c) + (mpfr_exp_t) (m - 3);
1774a238c70SJohn Marino if (!mpfr_can_round (c, err, MPFR_RNDN, MPFR_RNDZ,
1784a238c70SJohn Marino MPFR_PREC (z) + (rnd_mode == MPFR_RNDN)))
1794a238c70SJohn Marino goto next_step;
1804a238c70SJohn Marino
1814a238c70SJohn Marino /* we can't set z now, because in case z = x, and the mpfr_can_round()
1824a238c70SJohn Marino call below fails, we will have clobbered the input */
1834a238c70SJohn Marino mpfr_set_prec (xr, MPFR_PREC(c));
1844a238c70SJohn Marino mpfr_swap (xr, c); /* save the approximation of the cosine in xr */
1854a238c70SJohn Marino mpfr_sqr (c, xr, MPFR_RNDU); /* the absolute error is bounded by
1864a238c70SJohn Marino 2^(5-m) if reduce=1, and by 2^(2-m)
1874a238c70SJohn Marino otherwise */
1884a238c70SJohn Marino mpfr_ui_sub (c, 1, c, MPFR_RNDN); /* error bounded by 2^(6-m) if reduce
1894a238c70SJohn Marino is 1, and 2^(3-m) otherwise */
1904a238c70SJohn Marino mpfr_sqrt (c, c, MPFR_RNDN); /* the absolute error is bounded by
1914a238c70SJohn Marino 2^(6-m-Exp(c)) if reduce=1, and
1924a238c70SJohn Marino 2^(3-m-Exp(c)) otherwise */
1934a238c70SJohn Marino err = 3 + 3 * reduce - MPFR_GET_EXP (c);
1944a238c70SJohn Marino if (neg)
1954a238c70SJohn Marino MPFR_CHANGE_SIGN (c);
1964a238c70SJohn Marino
1974a238c70SJohn Marino /* the absolute error on c is at most 2^(err-m), which we must put
1984a238c70SJohn Marino in the form 2^(EXP(c)-err). */
1994a238c70SJohn Marino err = MPFR_GET_EXP (c) + (mpfr_exp_t) m - err;
2004a238c70SJohn Marino if (mpfr_can_round (c, err, MPFR_RNDN, MPFR_RNDZ,
2014a238c70SJohn Marino MPFR_PREC (y) + (rnd_mode == MPFR_RNDN)))
2024a238c70SJohn Marino break;
2034a238c70SJohn Marino /* check for huge cancellation */
2044a238c70SJohn Marino if (err < (mpfr_exp_t) MPFR_PREC (y))
2054a238c70SJohn Marino m += MPFR_PREC (y) - err;
2064a238c70SJohn Marino /* Check if near 1 */
2074a238c70SJohn Marino if (MPFR_GET_EXP (c) == 1
2084a238c70SJohn Marino && MPFR_MANT (c)[MPFR_LIMB_SIZE (c)-1] == MPFR_LIMB_HIGHBIT)
2094a238c70SJohn Marino m += m;
2104a238c70SJohn Marino
2114a238c70SJohn Marino next_step:
2124a238c70SJohn Marino MPFR_ZIV_NEXT (loop, m);
2134a238c70SJohn Marino mpfr_set_prec (c, m);
2144a238c70SJohn Marino }
2154a238c70SJohn Marino MPFR_ZIV_FREE (loop);
2164a238c70SJohn Marino
2174a238c70SJohn Marino inexy = mpfr_set (y, c, rnd_mode);
2184a238c70SJohn Marino inexz = mpfr_set (z, xr, rnd_mode);
2194a238c70SJohn Marino
2204a238c70SJohn Marino mpfr_clear (c);
2214a238c70SJohn Marino mpfr_clear (xr);
2224a238c70SJohn Marino
2234a238c70SJohn Marino end:
2244a238c70SJohn Marino MPFR_SAVE_EXPO_FREE (expo);
2254a238c70SJohn Marino /* FIXME: add a test for bug before revision 7355 */
2264a238c70SJohn Marino inexy = mpfr_check_range (y, inexy, rnd_mode);
2274a238c70SJohn Marino inexz = mpfr_check_range (z, inexz, rnd_mode);
2284a238c70SJohn Marino MPFR_RET (INEX(inexy,inexz));
2294a238c70SJohn Marino }
2304a238c70SJohn Marino
2314a238c70SJohn Marino /*************** asymptotically fast implementation below ********************/
2324a238c70SJohn Marino
2334a238c70SJohn Marino /* truncate Q from R to at most prec bits.
2344a238c70SJohn Marino Return the number of truncated bits.
2354a238c70SJohn Marino */
2364a238c70SJohn Marino static mpfr_prec_t
reduce(mpz_t Q,mpz_srcptr R,mpfr_prec_t prec)2374a238c70SJohn Marino reduce (mpz_t Q, mpz_srcptr R, mpfr_prec_t prec)
2384a238c70SJohn Marino {
2394a238c70SJohn Marino mpfr_prec_t l = mpz_sizeinbase (R, 2);
2404a238c70SJohn Marino
2414a238c70SJohn Marino l = (l > prec) ? l - prec : 0;
2424a238c70SJohn Marino mpz_fdiv_q_2exp (Q, R, l);
2434a238c70SJohn Marino return l;
2444a238c70SJohn Marino }
2454a238c70SJohn Marino
2464a238c70SJohn Marino /* truncate S and C so that the smaller has prec bits.
2474a238c70SJohn Marino Return the number of truncated bits.
2484a238c70SJohn Marino */
2494a238c70SJohn Marino static unsigned long
reduce2(mpz_t S,mpz_t C,mpfr_prec_t prec)2504a238c70SJohn Marino reduce2 (mpz_t S, mpz_t C, mpfr_prec_t prec)
2514a238c70SJohn Marino {
2524a238c70SJohn Marino unsigned long ls = mpz_sizeinbase (S, 2);
2534a238c70SJohn Marino unsigned long lc = mpz_sizeinbase (C, 2);
2544a238c70SJohn Marino unsigned long l;
2554a238c70SJohn Marino
2564a238c70SJohn Marino l = (ls < lc) ? ls : lc; /* smaller length */
2574a238c70SJohn Marino l = (l > prec) ? l - prec : 0;
2584a238c70SJohn Marino mpz_fdiv_q_2exp (S, S, l);
2594a238c70SJohn Marino mpz_fdiv_q_2exp (C, C, l);
2604a238c70SJohn Marino return l;
2614a238c70SJohn Marino }
2624a238c70SJohn Marino
2634a238c70SJohn Marino /* return in S0/Q0 a rational approximation of sin(X) with absolute error
2644a238c70SJohn Marino bounded by 9*2^(-prec), where 0 <= X=p/2^r <= 1/2,
2654a238c70SJohn Marino and in C0/Q0 a rational approximation of cos(X), with relative error
2664a238c70SJohn Marino bounded by 9*2^(-prec) (and also absolute error, since
2674a238c70SJohn Marino |cos(X)| <= 1).
2684a238c70SJohn Marino We have sin(X)/X = sum((-1)^i*(p/2^r)^i/(2i+1)!, i=0..infinity).
2694a238c70SJohn Marino We use the following binary splitting formula:
2704a238c70SJohn Marino P(a,b) = (-p)^(b-a)
2714a238c70SJohn Marino Q(a,b) = (2a)*(2a+1)*2^r if a+1=b [except Q(0,1)=1], Q(a,c)*Q(c,b) otherwise
2724a238c70SJohn Marino T(a,b) = 1 if a+1=b, Q(c,b)*T(a,c)+P(a,c)*T(c,b) otherwise.
2734a238c70SJohn Marino
2744a238c70SJohn Marino Since we use P(a,b) for b-a=2^k only, we compute only p^(2^k).
2754a238c70SJohn Marino We do not store the factor 2^r in Q().
2764a238c70SJohn Marino
2774a238c70SJohn Marino Then sin(X)/X ~ T(0,i)/Q(0,i) for i so that (p/2^r)^i/i! is small enough.
2784a238c70SJohn Marino
2794a238c70SJohn Marino Return l such that Q0 has to be multiplied by 2^l.
2804a238c70SJohn Marino
2814a238c70SJohn Marino Assumes prec >= 10.
2824a238c70SJohn Marino */
2834a238c70SJohn Marino static unsigned long
sin_bs_aux(mpz_t Q0,mpz_t S0,mpz_t C0,mpz_srcptr p,mpfr_prec_t r,mpfr_prec_t prec)2844a238c70SJohn Marino sin_bs_aux (mpz_t Q0, mpz_t S0, mpz_t C0, mpz_srcptr p, mpfr_prec_t r,
2854a238c70SJohn Marino mpfr_prec_t prec)
2864a238c70SJohn Marino {
2874a238c70SJohn Marino mpz_t T[GMP_NUMB_BITS], Q[GMP_NUMB_BITS], ptoj[GMP_NUMB_BITS], pp;
2884a238c70SJohn Marino mpfr_prec_t log2_nb_terms[GMP_NUMB_BITS], mult[GMP_NUMB_BITS];
2894a238c70SJohn Marino mpfr_prec_t accu[GMP_NUMB_BITS], size_ptoj[GMP_NUMB_BITS];
2904a238c70SJohn Marino mpfr_prec_t prec_i_have, r0 = r;
2914a238c70SJohn Marino unsigned long alloc, i, j, k;
2924a238c70SJohn Marino mpfr_prec_t l;
2934a238c70SJohn Marino
2944a238c70SJohn Marino if (MPFR_UNLIKELY(mpz_cmp_ui (p, 0) == 0)) /* sin(x)/x -> 1 */
2954a238c70SJohn Marino {
2964a238c70SJohn Marino mpz_set_ui (Q0, 1);
2974a238c70SJohn Marino mpz_set_ui (S0, 1);
2984a238c70SJohn Marino mpz_set_ui (C0, 1);
2994a238c70SJohn Marino return 0;
3004a238c70SJohn Marino }
3014a238c70SJohn Marino
3024a238c70SJohn Marino /* check that X=p/2^r <= 1/2 */
3034a238c70SJohn Marino MPFR_ASSERTN(mpz_sizeinbase (p, 2) - (mpfr_exp_t) r <= -1);
3044a238c70SJohn Marino
3054a238c70SJohn Marino mpz_init (pp);
3064a238c70SJohn Marino
3074a238c70SJohn Marino /* normalize p (non-zero here) */
3084a238c70SJohn Marino l = mpz_scan1 (p, 0);
3094a238c70SJohn Marino mpz_fdiv_q_2exp (pp, p, l); /* p = pp * 2^l */
3104a238c70SJohn Marino mpz_mul (pp, pp, pp);
3114a238c70SJohn Marino r = 2 * (r - l); /* x^2 = (p/2^r0)^2 = pp / 2^r */
3124a238c70SJohn Marino
3134a238c70SJohn Marino /* now p is odd */
3144a238c70SJohn Marino alloc = 2;
3154a238c70SJohn Marino mpz_init_set_ui (T[0], 6);
3164a238c70SJohn Marino mpz_init_set_ui (Q[0], 6);
3174a238c70SJohn Marino mpz_init_set (ptoj[0], pp); /* ptoj[i] = pp^(2^i) */
3184a238c70SJohn Marino mpz_init (T[1]);
3194a238c70SJohn Marino mpz_init (Q[1]);
3204a238c70SJohn Marino mpz_init (ptoj[1]);
3214a238c70SJohn Marino mpz_mul (ptoj[1], pp, pp); /* ptoj[1] = pp^2 */
3224a238c70SJohn Marino size_ptoj[1] = mpz_sizeinbase (ptoj[1], 2);
3234a238c70SJohn Marino
3244a238c70SJohn Marino mpz_mul_2exp (T[0], T[0], r);
3254a238c70SJohn Marino mpz_sub (T[0], T[0], pp); /* 6*2^r - pp = 6*2^r*(1 - x^2/6) */
3264a238c70SJohn Marino log2_nb_terms[0] = 1;
3274a238c70SJohn Marino
3284a238c70SJohn Marino /* already take into account the factor x=p/2^r in sin(x) = x * (...) */
3294a238c70SJohn Marino mult[0] = r - mpz_sizeinbase (pp, 2) + r0 - mpz_sizeinbase (p, 2);
3304a238c70SJohn Marino /* we have x^3 < 1/2^mult[0] */
3314a238c70SJohn Marino
3324a238c70SJohn Marino for (i = 2, k = 0, prec_i_have = mult[0]; prec_i_have < prec; i += 2)
3334a238c70SJohn Marino {
3344a238c70SJohn Marino /* i is even here */
3354a238c70SJohn Marino /* invariant: Q[0]*Q[1]*...*Q[k] equals (2i-1)!,
3364a238c70SJohn Marino we have already summed terms of index < i
3374a238c70SJohn Marino in S[0]/Q[0], ..., S[k]/Q[k] */
3384a238c70SJohn Marino k ++;
3394a238c70SJohn Marino if (k + 1 >= alloc) /* necessarily k + 1 = alloc */
3404a238c70SJohn Marino {
3414a238c70SJohn Marino alloc ++;
3424a238c70SJohn Marino mpz_init (T[k+1]);
3434a238c70SJohn Marino mpz_init (Q[k+1]);
3444a238c70SJohn Marino mpz_init (ptoj[k+1]);
3454a238c70SJohn Marino mpz_mul (ptoj[k+1], ptoj[k], ptoj[k]); /* pp^(2^(k+1)) */
3464a238c70SJohn Marino size_ptoj[k+1] = mpz_sizeinbase (ptoj[k+1], 2);
3474a238c70SJohn Marino }
3484a238c70SJohn Marino /* for i even, we have Q[k] = (2*i)*(2*i+1), T[k] = 1,
3494a238c70SJohn Marino then Q[k+1] = (2*i+2)*(2*i+3), T[k+1] = 1,
3504a238c70SJohn Marino which reduces to T[k] = (2*i+2)*(2*i+3)*2^r-pp,
3514a238c70SJohn Marino Q[k] = (2*i)*(2*i+1)*(2*i+2)*(2*i+3). */
3524a238c70SJohn Marino log2_nb_terms[k] = 1;
3534a238c70SJohn Marino mpz_set_ui (Q[k], (2 * i + 2) * (2 * i + 3));
3544a238c70SJohn Marino mpz_mul_2exp (T[k], Q[k], r);
3554a238c70SJohn Marino mpz_sub (T[k], T[k], pp);
3564a238c70SJohn Marino mpz_mul_ui (Q[k], Q[k], (2 * i) * (2 * i + 1));
3574a238c70SJohn Marino /* the next term of the series is divided by Q[k] and multiplied
3584a238c70SJohn Marino by pp^2/2^(2r), thus the mult. factor < 1/2^mult[k] */
3594a238c70SJohn Marino mult[k] = mpz_sizeinbase (Q[k], 2) + 2 * r - size_ptoj[1] - 1;
3604a238c70SJohn Marino /* the absolute contribution of the next term is 1/2^accu[k] */
3614a238c70SJohn Marino accu[k] = (k == 0) ? mult[k] : mult[k] + accu[k-1];
3624a238c70SJohn Marino prec_i_have = accu[k]; /* the current term is < 1/2^accu[k] */
3634a238c70SJohn Marino j = (i + 2) / 2;
3644a238c70SJohn Marino l = 1;
3654a238c70SJohn Marino while ((j & 1) == 0) /* combine and reduce */
3664a238c70SJohn Marino {
3674a238c70SJohn Marino mpz_mul (T[k], T[k], ptoj[l]);
3684a238c70SJohn Marino mpz_mul (T[k-1], T[k-1], Q[k]);
3694a238c70SJohn Marino mpz_mul_2exp (T[k-1], T[k-1], r << l);
3704a238c70SJohn Marino mpz_add (T[k-1], T[k-1], T[k]);
3714a238c70SJohn Marino mpz_mul (Q[k-1], Q[k-1], Q[k]);
3724a238c70SJohn Marino log2_nb_terms[k-1] ++; /* number of terms in S[k-1]
3734a238c70SJohn Marino is a power of 2 by construction */
3744a238c70SJohn Marino prec_i_have = mpz_sizeinbase (Q[k], 2);
3754a238c70SJohn Marino mult[k-1] += prec_i_have + (r << l) - size_ptoj[l] - 1;
3764a238c70SJohn Marino accu[k-1] = (k == 1) ? mult[k-1] : mult[k-1] + accu[k-2];
3774a238c70SJohn Marino prec_i_have = accu[k-1];
3784a238c70SJohn Marino l ++;
3794a238c70SJohn Marino j >>= 1;
3804a238c70SJohn Marino k --;
3814a238c70SJohn Marino }
3824a238c70SJohn Marino }
3834a238c70SJohn Marino
3844a238c70SJohn Marino /* accumulate all products in T[0] and Q[0]. Warning: contrary to above,
3854a238c70SJohn Marino here we do not have log2_nb_terms[k-1] = log2_nb_terms[k]+1. */
3864a238c70SJohn Marino l = 0; /* number of accumulated terms in the right part T[k]/Q[k] */
3874a238c70SJohn Marino while (k > 0)
3884a238c70SJohn Marino {
3894a238c70SJohn Marino j = log2_nb_terms[k-1];
3904a238c70SJohn Marino mpz_mul (T[k], T[k], ptoj[j]);
3914a238c70SJohn Marino mpz_mul (T[k-1], T[k-1], Q[k]);
3924a238c70SJohn Marino l += 1 << log2_nb_terms[k];
3934a238c70SJohn Marino mpz_mul_2exp (T[k-1], T[k-1], r * l);
3944a238c70SJohn Marino mpz_add (T[k-1], T[k-1], T[k]);
3954a238c70SJohn Marino mpz_mul (Q[k-1], Q[k-1], Q[k]);
3964a238c70SJohn Marino k--;
3974a238c70SJohn Marino }
3984a238c70SJohn Marino
3994a238c70SJohn Marino l = r0 + r * (i - 1); /* implicit multiplier 2^r for Q0 */
4004a238c70SJohn Marino /* at this point T[0]/(2^l*Q[0]) is an approximation of sin(x) where the 1st
4014a238c70SJohn Marino neglected term has contribution < 1/2^prec, thus since the series has
4024a238c70SJohn Marino alternate signs, the error is < 1/2^prec */
4034a238c70SJohn Marino
4044a238c70SJohn Marino /* we truncate Q0 to prec bits: the relative error is at most 2^(1-prec),
4054a238c70SJohn Marino which means that Q0 = Q[0] * (1+theta) with |theta| <= 2^(1-prec)
4064a238c70SJohn Marino [up to a power of two] */
4074a238c70SJohn Marino l += reduce (Q0, Q[0], prec);
4084a238c70SJohn Marino l -= reduce (T[0], T[0], prec);
4094a238c70SJohn Marino /* multiply by x = p/2^l */
4104a238c70SJohn Marino mpz_mul (S0, T[0], p);
4114a238c70SJohn Marino l -= reduce (S0, S0, prec); /* S0 = T[0] * (1 + theta)^2 up to power of 2 */
4124a238c70SJohn Marino /* sin(X) ~ S0/Q0*(1 + theta)^3 + err with |theta| <= 2^(1-prec) and
4134a238c70SJohn Marino |err| <= 2^(-prec), thus since |S0/Q0| <= 1:
4144a238c70SJohn Marino |sin(X) - S0/Q0| <= 4*|theta*S0/Q0| + |err| <= 9*2^(-prec) */
4154a238c70SJohn Marino
4164a238c70SJohn Marino mpz_clear (pp);
4174a238c70SJohn Marino for (j = 0; j < alloc; j ++)
4184a238c70SJohn Marino {
4194a238c70SJohn Marino mpz_clear (T[j]);
4204a238c70SJohn Marino mpz_clear (Q[j]);
4214a238c70SJohn Marino mpz_clear (ptoj[j]);
4224a238c70SJohn Marino }
4234a238c70SJohn Marino
4244a238c70SJohn Marino /* compute cos(X) from sin(X): sqrt(1-(S/Q)^2) = sqrt(Q^2-S^2)/Q
4254a238c70SJohn Marino = sqrt(Q0^2*2^(2l)-S0^2)/Q0.
4264a238c70SJohn Marino Write S/Q = sin(X) + eps with |eps| <= 9*2^(-prec),
4274a238c70SJohn Marino then sqrt(Q^2-S^2) = sqrt(Q^2-Q^2*(sin(X)+eps)^2)
4284a238c70SJohn Marino = sqrt(Q^2*cos(X)^2-Q^2*(2*sin(X)*eps+eps^2))
4294a238c70SJohn Marino = sqrt(Q^2*cos(X)^2-Q^2*eps1) with |eps1|<=9*2^(-prec)
4304a238c70SJohn Marino [using X<=1/2 and eps<=9*2^(-prec) and prec>=10]
4314a238c70SJohn Marino
4324a238c70SJohn Marino Since we truncate the square root, we get:
4334a238c70SJohn Marino sqrt(Q^2*cos(X)^2-Q^2*eps1)+eps2 with |eps2|<1
4344a238c70SJohn Marino = Q*sqrt(cos(X)^2-eps1)+eps2
4354a238c70SJohn Marino = Q*cos(X)*(1+eps3)+eps2 with |eps3| <= 6*2^(-prec)
4364a238c70SJohn Marino = Q*cos(X)*(1+eps3+eps2/(Q*cos(X)))
4374a238c70SJohn Marino = Q*cos(X)*(1+eps4) with |eps4| <= 9*2^(-prec)
4384a238c70SJohn Marino since |Q| >= 2^(prec-1) */
4394a238c70SJohn Marino /* we assume that Q0*2^l >= 2^(prec-1) */
4404a238c70SJohn Marino MPFR_ASSERTN(l + mpz_sizeinbase (Q0, 2) >= prec);
4414a238c70SJohn Marino mpz_mul (C0, Q0, Q0);
4424a238c70SJohn Marino mpz_mul_2exp (C0, C0, 2 * l);
4434a238c70SJohn Marino mpz_submul (C0, S0, S0);
4444a238c70SJohn Marino mpz_sqrt (C0, C0);
4454a238c70SJohn Marino
4464a238c70SJohn Marino return l;
4474a238c70SJohn Marino }
4484a238c70SJohn Marino
4494a238c70SJohn Marino /* Put in s and c approximations of sin(x) and cos(x) respectively.
4504a238c70SJohn Marino Assumes 0 < x < Pi/4 and PREC(s) = PREC(c) >= 10.
4514a238c70SJohn Marino Return err such that the relative error is bounded by 2^err ulps.
4524a238c70SJohn Marino */
4534a238c70SJohn Marino static int
sincos_aux(mpfr_t s,mpfr_t c,mpfr_srcptr x,mpfr_rnd_t rnd_mode)4544a238c70SJohn Marino sincos_aux (mpfr_t s, mpfr_t c, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
4554a238c70SJohn Marino {
4564a238c70SJohn Marino mpfr_prec_t prec_s, sh;
4574a238c70SJohn Marino mpz_t Q, S, C, Q2, S2, C2, y;
4584a238c70SJohn Marino mpfr_t x2;
4594a238c70SJohn Marino unsigned long l, l2, j, err;
4604a238c70SJohn Marino
4614a238c70SJohn Marino MPFR_ASSERTD(MPFR_PREC(s) == MPFR_PREC(c));
4624a238c70SJohn Marino
4634a238c70SJohn Marino prec_s = MPFR_PREC(s);
4644a238c70SJohn Marino
4654a238c70SJohn Marino mpfr_init2 (x2, MPFR_PREC(x));
4664a238c70SJohn Marino mpz_init (Q);
4674a238c70SJohn Marino mpz_init (S);
4684a238c70SJohn Marino mpz_init (C);
4694a238c70SJohn Marino mpz_init (Q2);
4704a238c70SJohn Marino mpz_init (S2);
4714a238c70SJohn Marino mpz_init (C2);
4724a238c70SJohn Marino mpz_init (y);
4734a238c70SJohn Marino
4744a238c70SJohn Marino mpfr_set (x2, x, MPFR_RNDN); /* exact */
4754a238c70SJohn Marino mpz_set_ui (Q, 1);
4764a238c70SJohn Marino l = 0;
4774a238c70SJohn Marino mpz_set_ui (S, 0); /* sin(0) = S/(2^l*Q), exact */
4784a238c70SJohn Marino mpz_set_ui (C, 1); /* cos(0) = C/(2^l*Q), exact */
4794a238c70SJohn Marino
4804a238c70SJohn Marino /* Invariant: x = X + x2/2^(sh-1), where the part X was already treated,
4814a238c70SJohn Marino S/(2^l*Q) ~ sin(X), C/(2^l*Q) ~ cos(X), and x2/2^(sh-1) < Pi/4.
4824a238c70SJohn Marino 'sh-1' is the number of already shifted bits in x2.
4834a238c70SJohn Marino */
4844a238c70SJohn Marino
4854a238c70SJohn Marino for (sh = 1, j = 0; mpfr_cmp_ui (x2, 0) != 0 && sh <= prec_s; sh <<= 1, j++)
4864a238c70SJohn Marino {
4874a238c70SJohn Marino if (sh > prec_s / 2) /* sin(x) = x + O(x^3), cos(x) = 1 + O(x^2) */
4884a238c70SJohn Marino {
4894a238c70SJohn Marino l2 = -mpfr_get_z_2exp (S2, x2); /* S2/2^l2 = x2 */
4904a238c70SJohn Marino l2 += sh - 1;
4914a238c70SJohn Marino mpz_set_ui (Q2, 1);
4924a238c70SJohn Marino mpz_set_ui (C2, 1);
4934a238c70SJohn Marino mpz_mul_2exp (C2, C2, l2);
4944a238c70SJohn Marino mpfr_set_ui (x2, 0, MPFR_RNDN);
4954a238c70SJohn Marino }
4964a238c70SJohn Marino else
4974a238c70SJohn Marino {
4984a238c70SJohn Marino /* y <- trunc(x2 * 2^sh) = trunc(x * 2^(2*sh-1)) */
4994a238c70SJohn Marino mpfr_mul_2exp (x2, x2, sh, MPFR_RNDN); /* exact */
5004a238c70SJohn Marino mpfr_get_z (y, x2, MPFR_RNDZ); /* round towards zero: now
5014a238c70SJohn Marino 0 <= x2 < 2^sh, thus
5024a238c70SJohn Marino 0 <= x2/2^(sh-1) < 2^(1-sh) */
5034a238c70SJohn Marino if (mpz_cmp_ui (y, 0) == 0)
5044a238c70SJohn Marino continue;
5054a238c70SJohn Marino mpfr_sub_z (x2, x2, y, MPFR_RNDN); /* should be exact */
5064a238c70SJohn Marino l2 = sin_bs_aux (Q2, S2, C2, y, 2 * sh - 1, prec_s);
5074a238c70SJohn Marino /* we now have |S2/Q2/2^l2 - sin(X)| <= 9*2^(prec_s)
5084a238c70SJohn Marino and |C2/Q2/2^l2 - cos(X)| <= 6*2^(prec_s), with X=y/2^(2sh-1) */
5094a238c70SJohn Marino }
5104a238c70SJohn Marino if (sh == 1) /* S=0, C=1 */
5114a238c70SJohn Marino {
5124a238c70SJohn Marino l = l2;
5134a238c70SJohn Marino mpz_swap (Q, Q2);
5144a238c70SJohn Marino mpz_swap (S, S2);
5154a238c70SJohn Marino mpz_swap (C, C2);
5164a238c70SJohn Marino }
5174a238c70SJohn Marino else
5184a238c70SJohn Marino {
5194a238c70SJohn Marino /* s <- s*c2+c*s2, c <- c*c2-s*s2, using Karatsuba:
5204a238c70SJohn Marino a = s+c, b = s2+c2, t = a*b, d = s*s2, e = c*c2,
5214a238c70SJohn Marino s <- t - d - e, c <- e - d */
5224a238c70SJohn Marino mpz_add (y, S, C); /* a */
5234a238c70SJohn Marino mpz_mul (C, C, C2); /* e */
5244a238c70SJohn Marino mpz_add (C2, C2, S2); /* b */
5254a238c70SJohn Marino mpz_mul (S2, S, S2); /* d */
5264a238c70SJohn Marino mpz_mul (y, y, C2); /* a*b */
5274a238c70SJohn Marino mpz_sub (S, y, S2); /* t - d */
5284a238c70SJohn Marino mpz_sub (S, S, C); /* t - d - e */
5294a238c70SJohn Marino mpz_sub (C, C, S2); /* e - d */
5304a238c70SJohn Marino mpz_mul (Q, Q, Q2);
5314a238c70SJohn Marino /* after j loops, the error is <= (11j-2)*2^(prec_s) */
5324a238c70SJohn Marino l += l2;
5334a238c70SJohn Marino /* reduce Q to prec_s bits */
5344a238c70SJohn Marino l += reduce (Q, Q, prec_s);
5354a238c70SJohn Marino /* reduce S,C to prec_s bits, error <= 11*j*2^(prec_s) */
5364a238c70SJohn Marino l -= reduce2 (S, C, prec_s);
5374a238c70SJohn Marino }
5384a238c70SJohn Marino }
5394a238c70SJohn Marino
5404a238c70SJohn Marino j = 11 * j;
5414a238c70SJohn Marino for (err = 0; j > 1; j = (j + 1) / 2, err ++);
5424a238c70SJohn Marino
5434a238c70SJohn Marino mpfr_set_z (s, S, MPFR_RNDN);
5444a238c70SJohn Marino mpfr_div_z (s, s, Q, MPFR_RNDN);
5454a238c70SJohn Marino mpfr_div_2exp (s, s, l, MPFR_RNDN);
5464a238c70SJohn Marino
5474a238c70SJohn Marino mpfr_set_z (c, C, MPFR_RNDN);
5484a238c70SJohn Marino mpfr_div_z (c, c, Q, MPFR_RNDN);
5494a238c70SJohn Marino mpfr_div_2exp (c, c, l, MPFR_RNDN);
5504a238c70SJohn Marino
5514a238c70SJohn Marino mpz_clear (Q);
5524a238c70SJohn Marino mpz_clear (S);
5534a238c70SJohn Marino mpz_clear (C);
5544a238c70SJohn Marino mpz_clear (Q2);
5554a238c70SJohn Marino mpz_clear (S2);
5564a238c70SJohn Marino mpz_clear (C2);
5574a238c70SJohn Marino mpz_clear (y);
5584a238c70SJohn Marino mpfr_clear (x2);
5594a238c70SJohn Marino return err;
5604a238c70SJohn Marino }
5614a238c70SJohn Marino
5624a238c70SJohn Marino /* Assumes x is neither NaN, +/-Inf, nor +/- 0.
5634a238c70SJohn Marino One of s and c might be NULL, in which case the corresponding value is
5644a238c70SJohn Marino not computed.
5654a238c70SJohn Marino Assumes s differs from c.
5664a238c70SJohn Marino */
5674a238c70SJohn Marino int
mpfr_sincos_fast(mpfr_t s,mpfr_t c,mpfr_srcptr x,mpfr_rnd_t rnd)5684a238c70SJohn Marino mpfr_sincos_fast (mpfr_t s, mpfr_t c, mpfr_srcptr x, mpfr_rnd_t rnd)
5694a238c70SJohn Marino {
5704a238c70SJohn Marino int inexs, inexc;
5714a238c70SJohn Marino mpfr_t x_red, ts, tc;
5724a238c70SJohn Marino mpfr_prec_t w;
5734a238c70SJohn Marino mpfr_exp_t err, errs, errc;
5744a238c70SJohn Marino MPFR_ZIV_DECL (loop);
5754a238c70SJohn Marino
5764a238c70SJohn Marino MPFR_ASSERTN(s != c);
5774a238c70SJohn Marino if (s == NULL)
5784a238c70SJohn Marino w = MPFR_PREC(c);
5794a238c70SJohn Marino else if (c == NULL)
5804a238c70SJohn Marino w = MPFR_PREC(s);
5814a238c70SJohn Marino else
5824a238c70SJohn Marino w = MPFR_PREC(s) >= MPFR_PREC(c) ? MPFR_PREC(s) : MPFR_PREC(c);
5834a238c70SJohn Marino w += MPFR_INT_CEIL_LOG2(w) + 9; /* ensures w >= 10 (needed by sincos_aux) */
5844a238c70SJohn Marino mpfr_init2 (ts, w);
5854a238c70SJohn Marino mpfr_init2 (tc, w);
5864a238c70SJohn Marino
5874a238c70SJohn Marino MPFR_ZIV_INIT (loop, w);
5884a238c70SJohn Marino for (;;)
5894a238c70SJohn Marino {
5904a238c70SJohn Marino /* if 0 < x <= Pi/4, we can call sincos_aux directly */
5914a238c70SJohn Marino if (MPFR_IS_POS(x) && mpfr_cmp_ui_2exp (x, 1686629713, -31) <= 0)
5924a238c70SJohn Marino {
5934a238c70SJohn Marino err = sincos_aux (ts, tc, x, MPFR_RNDN);
5944a238c70SJohn Marino }
5954a238c70SJohn Marino /* if -Pi/4 <= x < 0, use sin(-x)=-sin(x) */
5964a238c70SJohn Marino else if (MPFR_IS_NEG(x) && mpfr_cmp_si_2exp (x, -1686629713, -31) >= 0)
5974a238c70SJohn Marino {
5984a238c70SJohn Marino mpfr_init2 (x_red, MPFR_PREC(x));
5994a238c70SJohn Marino mpfr_neg (x_red, x, rnd); /* exact */
6004a238c70SJohn Marino err = sincos_aux (ts, tc, x_red, MPFR_RNDN);
6014a238c70SJohn Marino mpfr_neg (ts, ts, MPFR_RNDN);
6024a238c70SJohn Marino mpfr_clear (x_red);
6034a238c70SJohn Marino }
6044a238c70SJohn Marino else /* argument reduction is needed */
6054a238c70SJohn Marino {
6064a238c70SJohn Marino long q;
6074a238c70SJohn Marino mpfr_t pi;
6084a238c70SJohn Marino int neg = 0;
6094a238c70SJohn Marino
6104a238c70SJohn Marino mpfr_init2 (x_red, w);
6114a238c70SJohn Marino mpfr_init2 (pi, (MPFR_EXP(x) > 0) ? w + MPFR_EXP(x) : w);
6124a238c70SJohn Marino mpfr_const_pi (pi, MPFR_RNDN);
6134a238c70SJohn Marino mpfr_div_2exp (pi, pi, 1, MPFR_RNDN); /* Pi/2 */
6144a238c70SJohn Marino mpfr_remquo (x_red, &q, x, pi, MPFR_RNDN);
6154a238c70SJohn Marino /* x = q * (Pi/2 + eps1) + x_red + eps2,
6164a238c70SJohn Marino where |eps1| <= 1/2*ulp(Pi/2) = 2^(-w-MAX(0,EXP(x))),
6174a238c70SJohn Marino and eps2 <= 1/2*ulp(x_red) <= 1/2*ulp(Pi/2) = 2^(-w)
6184a238c70SJohn Marino Since |q| <= x/(Pi/2) <= |x|, we have
6194a238c70SJohn Marino q*|eps1| <= 2^(-w), thus
6204a238c70SJohn Marino |x - q * Pi/2 - x_red| <= 2^(1-w) */
6214a238c70SJohn Marino /* now -Pi/4 <= x_red <= Pi/4: if x_red < 0, consider -x_red */
6224a238c70SJohn Marino if (MPFR_IS_NEG(x_red))
6234a238c70SJohn Marino {
6244a238c70SJohn Marino mpfr_neg (x_red, x_red, MPFR_RNDN);
6254a238c70SJohn Marino neg = 1;
6264a238c70SJohn Marino }
6274a238c70SJohn Marino err = sincos_aux (ts, tc, x_red, MPFR_RNDN);
6284a238c70SJohn Marino err ++; /* to take into account the argument reduction */
6294a238c70SJohn Marino if (neg) /* sin(-x) = -sin(x), cos(-x) = cos(x) */
6304a238c70SJohn Marino mpfr_neg (ts, ts, MPFR_RNDN);
6314a238c70SJohn Marino if (q & 2) /* sin(x+Pi) = -sin(x), cos(x+Pi) = -cos(x) */
6324a238c70SJohn Marino {
6334a238c70SJohn Marino mpfr_neg (ts, ts, MPFR_RNDN);
6344a238c70SJohn Marino mpfr_neg (tc, tc, MPFR_RNDN);
6354a238c70SJohn Marino }
6364a238c70SJohn Marino if (q & 1) /* sin(x+Pi/2) = cos(x), cos(x+Pi/2) = -sin(x) */
6374a238c70SJohn Marino {
6384a238c70SJohn Marino mpfr_neg (ts, ts, MPFR_RNDN);
6394a238c70SJohn Marino mpfr_swap (ts, tc);
6404a238c70SJohn Marino }
6414a238c70SJohn Marino mpfr_clear (x_red);
6424a238c70SJohn Marino mpfr_clear (pi);
6434a238c70SJohn Marino }
6444a238c70SJohn Marino /* adjust errors with respect to absolute values */
6454a238c70SJohn Marino errs = err - MPFR_EXP(ts);
6464a238c70SJohn Marino errc = err - MPFR_EXP(tc);
6474a238c70SJohn Marino if ((s == NULL || MPFR_CAN_ROUND (ts, w - errs, MPFR_PREC(s), rnd)) &&
6484a238c70SJohn Marino (c == NULL || MPFR_CAN_ROUND (tc, w - errc, MPFR_PREC(c), rnd)))
6494a238c70SJohn Marino break;
6504a238c70SJohn Marino MPFR_ZIV_NEXT (loop, w);
6514a238c70SJohn Marino mpfr_set_prec (ts, w);
6524a238c70SJohn Marino mpfr_set_prec (tc, w);
6534a238c70SJohn Marino }
6544a238c70SJohn Marino MPFR_ZIV_FREE (loop);
6554a238c70SJohn Marino
6564a238c70SJohn Marino inexs = (s == NULL) ? 0 : mpfr_set (s, ts, rnd);
6574a238c70SJohn Marino inexc = (c == NULL) ? 0 : mpfr_set (c, tc, rnd);
6584a238c70SJohn Marino
6594a238c70SJohn Marino mpfr_clear (ts);
6604a238c70SJohn Marino mpfr_clear (tc);
6614a238c70SJohn Marino return INEX(inexs,inexc);
6624a238c70SJohn Marino }
663