Home | History | Annotate | Line # | Download | only in src
      1 /* mpc_sin_cos -- combined sine and cosine of a complex number.
      2 
      3 Copyright (C) 2010, 2011, 2012, 2020 INRIA
      4 
      5 This file is part of GNU MPC.
      6 
      7 GNU MPC is free software; you can redistribute it and/or modify it under
      8 the terms of the GNU Lesser General Public License as published by the
      9 Free Software Foundation; either version 3 of the License, or (at your
     10 option) any later version.
     11 
     12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
     13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
     14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
     15 more details.
     16 
     17 You should have received a copy of the GNU Lesser General Public License
     18 along with this program. If not, see http://www.gnu.org/licenses/ .
     19 */
     20 
     21 #include <stdio.h>
     22 #include "mpc-impl.h"
     23 
     24 static int
     25 mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
     26    mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
     27    /* assumes that op (that is, its real or imaginary part) is not finite */
     28 {
     29    int overlap;
     30    mpc_t op_loc;
     31 
     32    overlap = (rop_sin == op || rop_cos == op);
     33    if (overlap) {
     34       mpc_init3 (op_loc, MPC_PREC_RE (op), MPC_PREC_IM (op));
     35       mpc_set (op_loc, op, MPC_RNDNN);
     36    }
     37    else
     38       op_loc [0] = op [0];
     39 
     40    if (rop_sin != NULL) {
     41       if (mpfr_nan_p (mpc_realref (op_loc)) || mpfr_nan_p (mpc_imagref (op_loc))) {
     42          mpc_set (rop_sin, op_loc, rnd_sin);
     43          if (mpfr_nan_p (mpc_imagref (op_loc))) {
     44             /* sin(x +i*NaN) = NaN +i*NaN, except for x=0 */
     45             /* sin(-0 +i*NaN) = -0 +i*NaN */
     46             /* sin(+0 +i*NaN) = +0 +i*NaN */
     47             if (!mpfr_zero_p (mpc_realref (op_loc)))
     48                mpfr_set_nan (mpc_realref (rop_sin));
     49          }
     50          else /* op = NaN + i*y */
     51             if (!mpfr_inf_p (mpc_imagref (op_loc)) && !mpfr_zero_p (mpc_imagref (op_loc)))
     52             /* sin(NaN -i*Inf) = NaN -i*Inf */
     53             /* sin(NaN -i*0) = NaN -i*0 */
     54             /* sin(NaN +i*0) = NaN +i*0 */
     55             /* sin(NaN +i*Inf) = NaN +i*Inf */
     56             /* sin(NaN +i*y) = NaN +i*NaN, when 0<|y|<Inf */
     57             mpfr_set_nan (mpc_imagref (rop_sin));
     58       }
     59       else if (mpfr_inf_p (mpc_realref (op_loc))) {
     60          mpfr_set_nan (mpc_realref (rop_sin));
     61 
     62          if (!mpfr_inf_p (mpc_imagref (op_loc)) && !mpfr_zero_p (mpc_imagref (op_loc)))
     63             /* sin(+/-Inf +i*y) = NaN +i*NaN, when 0<|y|<Inf */
     64             mpfr_set_nan (mpc_imagref (rop_sin));
     65          else
     66             /* sin(+/-Inf -i*Inf) = NaN -i*Inf */
     67             /* sin(+/-Inf +i*Inf) = NaN +i*Inf */
     68             /* sin(+/-Inf -i*0) = NaN -i*0 */
     69             /* sin(+/-Inf +i*0) = NaN +i*0 */
     70             mpfr_set (mpc_imagref (rop_sin), mpc_imagref (op_loc), MPC_RND_IM (rnd_sin));
     71       }
     72       else if (mpfr_zero_p (mpc_realref (op_loc))) {
     73          /* sin(-0 -i*Inf) = -0 -i*Inf */
     74          /* sin(+0 -i*Inf) = +0 -i*Inf */
     75          /* sin(-0 +i*Inf) = -0 +i*Inf */
     76          /* sin(+0 +i*Inf) = +0 +i*Inf */
     77          mpc_set (rop_sin, op_loc, rnd_sin);
     78       }
     79       else {
     80          /* sin(x -i*Inf) = +Inf*(sin(x) -i*cos(x)) */
     81          /* sin(x +i*Inf) = +Inf*(sin(x) +i*cos(x)) */
     82          mpfr_t s, c;
     83          mpfr_init2 (s, 2);
     84          mpfr_init2 (c, 2);
     85          mpfr_sin_cos (s, c, mpc_realref (op_loc), MPFR_RNDZ);
     86          mpfr_set_inf (mpc_realref (rop_sin), MPFR_SIGN (s));
     87          mpfr_set_inf (mpc_imagref (rop_sin), MPFR_SIGN (c)*MPFR_SIGN (mpc_imagref (op_loc)));
     88          mpfr_clear (s);
     89          mpfr_clear (c);
     90       }
     91    }
     92 
     93    if (rop_cos != NULL) {
     94       if (mpfr_nan_p (mpc_realref (op_loc))) {
     95          /* cos(NaN + i * NaN) = NaN + i * NaN */
     96          /* cos(NaN - i * Inf) = +Inf + i * NaN */
     97          /* cos(NaN + i * Inf) = +Inf + i * NaN */
     98          /* cos(NaN - i * 0) = NaN - i * 0 */
     99          /* cos(NaN + i * 0) = NaN + i * 0 */
    100          /* cos(NaN + i * y) = NaN + i * NaN, when y != 0 */
    101          if (mpfr_inf_p (mpc_imagref (op_loc)))
    102             mpfr_set_inf (mpc_realref (rop_cos), +1);
    103          else
    104             mpfr_set_nan (mpc_realref (rop_cos));
    105 
    106          if (mpfr_zero_p (mpc_imagref (op_loc)))
    107             mpfr_set (mpc_imagref (rop_cos), mpc_imagref (op_loc), MPC_RND_IM (rnd_cos));
    108          else
    109             mpfr_set_nan (mpc_imagref (rop_cos));
    110       }
    111       else if (mpfr_nan_p (mpc_imagref (op_loc))) {
    112           /* cos(-Inf + i * NaN) = NaN + i * NaN */
    113           /* cos(+Inf + i * NaN) = NaN + i * NaN */
    114           /* cos(-0 + i * NaN) = NaN - i * 0 */
    115           /* cos(+0 + i * NaN) = NaN + i * 0 */
    116           /* cos(x + i * NaN) = NaN + i * NaN, when x != 0 */
    117          if (mpfr_zero_p (mpc_realref (op_loc)))
    118             mpfr_set (mpc_imagref (rop_cos), mpc_realref (op_loc), MPC_RND_IM (rnd_cos));
    119          else
    120             mpfr_set_nan (mpc_imagref (rop_cos));
    121 
    122          mpfr_set_nan (mpc_realref (rop_cos));
    123       }
    124       else if (mpfr_inf_p (mpc_realref (op_loc))) {
    125          /* cos(-Inf -i*Inf) = cos(+Inf +i*Inf) = -Inf +i*NaN */
    126          /* cos(-Inf +i*Inf) = cos(+Inf -i*Inf) = +Inf +i*NaN */
    127          /* cos(-Inf -i*0) = cos(+Inf +i*0) = NaN -i*0 */
    128          /* cos(-Inf +i*0) = cos(+Inf -i*0) = NaN +i*0 */
    129          /* cos(-Inf +i*y) = cos(+Inf +i*y) = NaN +i*NaN, when y != 0 */
    130 
    131          const int same_sign =
    132             mpfr_signbit (mpc_realref (op_loc)) == mpfr_signbit (mpc_imagref (op_loc));
    133 
    134          if (mpfr_inf_p (mpc_imagref (op_loc)))
    135             mpfr_set_inf (mpc_realref (rop_cos), (same_sign ? -1 : +1));
    136          else
    137             mpfr_set_nan (mpc_realref (rop_cos));
    138 
    139          if (mpfr_zero_p (mpc_imagref (op_loc)))
    140             mpfr_setsign (mpc_imagref (rop_cos), mpc_imagref (op_loc), same_sign,
    141                           MPC_RND_IM(rnd_cos));
    142          else
    143             mpfr_set_nan (mpc_imagref (rop_cos));
    144       }
    145       else if (mpfr_zero_p (mpc_realref (op_loc))) {
    146          /* cos(-0 -i*Inf) = cos(+0 +i*Inf) = +Inf -i*0 */
    147          /* cos(-0 +i*Inf) = cos(+0 -i*Inf) = +Inf +i*0 */
    148          const int same_sign =
    149             mpfr_signbit (mpc_realref (op_loc)) == mpfr_signbit (mpc_imagref (op_loc));
    150 
    151          mpfr_setsign (mpc_imagref (rop_cos), mpc_realref (op_loc), same_sign,
    152                        MPC_RND_IM (rnd_cos));
    153          mpfr_set_inf (mpc_realref (rop_cos), +1);
    154       }
    155       else {
    156          /* cos(x -i*Inf) = +Inf*cos(x) +i*Inf*sin(x), when x != 0 */
    157          /* cos(x +i*Inf) = +Inf*cos(x) -i*Inf*sin(x), when x != 0 */
    158          mpfr_t s, c;
    159          mpfr_init2 (c, 2);
    160          mpfr_init2 (s, 2);
    161          mpfr_sin_cos (s, c, mpc_realref (op_loc), MPFR_RNDN);
    162          mpfr_set_inf (mpc_realref (rop_cos), mpfr_sgn (c));
    163          mpfr_set_inf (mpc_imagref (rop_cos),
    164             (mpfr_sgn (mpc_imagref (op_loc)) == mpfr_sgn (s) ? -1 : +1));
    165          mpfr_clear (s);
    166          mpfr_clear (c);
    167       }
    168    }
    169 
    170    if (overlap)
    171       mpc_clear (op_loc);
    172 
    173    return MPC_INEX12 (MPC_INEX (0,0), MPC_INEX (0,0));
    174       /* everything is exact */
    175 }
    176 
    177 
    178 static int
    179 mpc_sin_cos_real (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
    180    mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
    181    /* assumes that op is real */
    182 {
    183    int inex_sin_re = 0, inex_cos_re = 0;
    184       /* Until further notice, assume computations exact; in particular,
    185          by definition, for not computed values.                         */
    186    mpfr_t s, c;
    187    int inex_s, inex_c;
    188    int sign_im = mpfr_signbit (mpc_imagref (op));
    189 
    190    /* sin(x +-0*i) = sin(x) +-0*i*sign(cos(x)) */
    191    /* cos(x +-i*0) = cos(x) -+i*0*sign(sin(x)) */
    192    if (rop_sin != 0)
    193       mpfr_init2 (s, MPC_PREC_RE (rop_sin));
    194    else
    195       mpfr_init2 (s, 2); /* We need only the sign. */
    196    if (rop_cos != NULL)
    197       mpfr_init2 (c, MPC_PREC_RE (rop_cos));
    198    else
    199       mpfr_init2 (c, 2);
    200    inex_s = mpfr_sin (s, mpc_realref (op), MPC_RND_RE (rnd_sin));
    201    inex_c = mpfr_cos (c, mpc_realref (op), MPC_RND_RE (rnd_cos));
    202       /* We cannot use mpfr_sin_cos since we may need two distinct rounding
    203          modes and the exact return values. If we need only the sign, an
    204          arbitrary rounding mode will work.                                 */
    205 
    206    if (rop_sin != NULL) {
    207       mpfr_set (mpc_realref (rop_sin), s, MPFR_RNDN); /* exact */
    208       inex_sin_re = inex_s;
    209       mpfr_set_zero (mpc_imagref (rop_sin),
    210          (     ( sign_im && !mpfr_signbit(c))
    211             || (!sign_im &&  mpfr_signbit(c)) ? -1 : 1));
    212    }
    213 
    214    if (rop_cos != NULL) {
    215       mpfr_set (mpc_realref (rop_cos), c, MPFR_RNDN); /* exact */
    216       inex_cos_re = inex_c;
    217       mpfr_set_zero (mpc_imagref (rop_cos),
    218          (     ( sign_im &&  mpfr_signbit(s))
    219             || (!sign_im && !mpfr_signbit(s)) ? -1 : 1));
    220    }
    221 
    222    mpfr_clear (s);
    223    mpfr_clear (c);
    224 
    225    return MPC_INEX12 (MPC_INEX (inex_sin_re, 0), MPC_INEX (inex_cos_re, 0));
    226 }
    227 
    228 
    229 static int
    230 mpc_sin_cos_imag (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
    231    mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
    232    /* assumes that op is purely imaginary, but not zero */
    233 {
    234    int inex_sin_im = 0, inex_cos_re = 0;
    235       /* assume exact if not computed */
    236    int overlap;
    237    mpc_t op_loc;
    238 
    239    overlap = (rop_sin == op || rop_cos == op);
    240    if (overlap) {
    241       mpc_init3 (op_loc, MPC_PREC_RE (op), MPC_PREC_IM (op));
    242       mpc_set (op_loc, op, MPC_RNDNN);
    243    }
    244    else
    245       op_loc [0] = op [0];
    246 
    247    if (rop_sin != NULL) {
    248       /* sin(+-O +i*y) = +-0 +i*sinh(y) */
    249       mpfr_set (mpc_realref(rop_sin), mpc_realref(op_loc), MPFR_RNDN);
    250       inex_sin_im = mpfr_sinh (mpc_imagref(rop_sin), mpc_imagref(op_loc), MPC_RND_IM(rnd_sin));
    251    }
    252 
    253    if (rop_cos != NULL) {
    254       /* cos(-0 - i * y) = cos(+0 + i * y) = cosh(y) - i * 0,
    255          cos(-0 + i * y) = cos(+0 - i * y) = cosh(y) + i * 0,
    256          where y > 0 */
    257       inex_cos_re = mpfr_cosh (mpc_realref (rop_cos), mpc_imagref (op_loc), MPC_RND_RE (rnd_cos));
    258 
    259       mpfr_set_ui (mpc_imagref (rop_cos), 0ul, MPC_RND_IM (rnd_cos));
    260       if (mpfr_signbit (mpc_realref (op_loc)) ==  mpfr_signbit (mpc_imagref (op_loc)))
    261          MPFR_CHANGE_SIGN (mpc_imagref (rop_cos));
    262    }
    263 
    264    if (overlap)
    265       mpc_clear (op_loc);
    266 
    267    return MPC_INEX12 (MPC_INEX (0, inex_sin_im), MPC_INEX (inex_cos_re, 0));
    268 }
    269 
    270 /* Fix an inexact overflow, when x is +Inf or -Inf:
    271    When rnd is towards zero, change x into the largest (in absolute value)
    272    floating-point number.
    273    Return the inexact flag. */
    274 int
    275 mpc_fix_inf (mpfr_t x, mpfr_rnd_t rnd)
    276 {
    277   MPC_ASSERT (mpfr_inf_p (x));
    278   if (!MPC_IS_LIKE_RNDZ(rnd, MPFR_SIGNBIT(x)))
    279     return mpfr_sgn (x);
    280   else
    281     {
    282       if (mpfr_sgn (x) < 0)
    283         mpfr_nextabove (x);
    284       else
    285         mpfr_nextbelow (x);
    286       return -mpfr_sgn (x);
    287     }
    288 }
    289 
    290 /* Fix an inexact underflow, when x is +0 or -0:
    291    When rnd is away from zero, change x into the closest floating-point number.
    292    Return the inexact flag. */
    293 int
    294 mpc_fix_zero (mpfr_t x, mpfr_rnd_t rnd)
    295 {
    296   if (!MPC_IS_LIKE_RNDA(rnd, MPFR_SIGNBIT(x)))
    297     return mpfr_signbit (x) == 0 ? -1 : 1;
    298   else
    299     {
    300       if (mpfr_signbit (x) == 0)
    301         {
    302           mpfr_nextabove (x);
    303           return 1;
    304         }
    305       else
    306         {
    307           mpfr_nextbelow (x);
    308           return -1;
    309         }
    310     }
    311 }
    312 
    313 int
    314 mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
    315    mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
    316    /* Feature not documented in the texinfo file: One of rop_sin or
    317       rop_cos may be NULL, in which case it is not computed, and the
    318       corresponding ternary inexact value is set to 0 (exact).       */
    319 {
    320    if (!mpc_fin_p (op))
    321       return mpc_sin_cos_nonfinite (rop_sin, rop_cos, op, rnd_sin, rnd_cos);
    322    else if (mpfr_zero_p (mpc_imagref (op)))
    323       return mpc_sin_cos_real (rop_sin, rop_cos, op, rnd_sin, rnd_cos);
    324    else if (mpfr_zero_p (mpc_realref (op)))
    325       return mpc_sin_cos_imag (rop_sin, rop_cos, op, rnd_sin, rnd_cos);
    326    else {
    327       /* let op = a + i*b, then sin(op) = sin(a)*cosh(b) + i*cos(a)*sinh(b)
    328                            and  cos(op) = cos(a)*cosh(b) - i*sin(a)*sinh(b).
    329 
    330          For Re(sin(op)) (and analogously, the other parts), we use the
    331          following algorithm, with rounding to nearest for all operations
    332          and working precision w:
    333 
    334          (1) x = o(sin(a))
    335          (2) y = o(cosh(b))
    336          (3) r = o(x*y)
    337          then the error on r is at most 4 ulps, since we can write
    338          r = sin(a)*cosh(b)*(1+t)^3 with |t| <= 2^(-w),
    339          thus for w >= 2, r = sin(a)*cosh(b)*(1+4*t) with |t| <= 2^(-w),
    340          thus the relative error is bounded by 4*2^(-w) <= 4*ulp(r).
    341       */
    342       mpfr_t s, c, sh, ch, sch, csh;
    343       mpfr_prec_t prec;
    344       int ok;
    345       int inex_re, inex_im, inex_sin, inex_cos, loop = 0;
    346       mpfr_exp_t saved_emin, saved_emax;
    347 
    348       saved_emin = mpfr_get_emin ();
    349       saved_emax = mpfr_get_emax ();
    350       mpfr_set_emin (mpfr_get_emin_min ());
    351       mpfr_set_emax (mpfr_get_emax_max ());
    352 
    353       prec = 2;
    354       if (rop_sin != NULL)
    355         {
    356           mp_prec_t er, ei;
    357           prec = MPC_MAX (prec, MPC_MAX_PREC (rop_sin));
    358           /* since the Taylor expansion of sin(x) at x=0 is x - x^3/6 + O(x^5),
    359              if x <= 2^(-p), then the second term/x is about 2^(-2p)/6, thus we
    360              need at least 2p+3 bits of precision. This is true only when x is
    361              exactly representable in the target precision. */
    362           if (MPC_MAX_PREC (op) <= prec)
    363             {
    364               er = mpfr_get_exp (mpc_realref (op));
    365               ei = mpfr_get_exp (mpc_imagref (op));
    366               /* consider the maximal exponent only */
    367               er = (er < ei) ? ei : er;
    368               if (er < 0)
    369                 if (prec < 2 * (mp_prec_t) (-er) + 3)
    370                   prec = 2 * (mp_prec_t) (-er) + 3;
    371             }
    372         }
    373       if (rop_cos != NULL)
    374          prec = MPC_MAX (prec, MPC_MAX_PREC (rop_cos));
    375 
    376       mpfr_init2 (s, 2);
    377       mpfr_init2 (c, 2);
    378       mpfr_init2 (sh, 2);
    379       mpfr_init2 (ch, 2);
    380       mpfr_init2 (sch, 2);
    381       mpfr_init2 (csh, 2);
    382 
    383       do {
    384          loop ++;
    385          prec += (loop <= 2) ? mpc_ceil_log2 (prec) + 5 : prec / 2;
    386 
    387          mpfr_set_prec (s, prec);
    388          mpfr_set_prec (c, prec);
    389          mpfr_set_prec (sh, prec);
    390          mpfr_set_prec (ch, prec);
    391          mpfr_set_prec (sch, prec);
    392          mpfr_set_prec (csh, prec);
    393 
    394          mpfr_sin_cos (s, c, mpc_realref(op), MPFR_RNDN);
    395          mpfr_sinh_cosh (sh, ch, mpc_imagref(op), MPFR_RNDN);
    396 
    397          ok = 1;
    398 
    399          if (rop_sin != NULL) {
    400             /* real part of sine */
    401             mpfr_mul (sch, s, ch, MPFR_RNDN);
    402             ok = (!mpfr_number_p (sch))
    403                   || mpfr_can_round (sch, prec - 2, MPFR_RNDN, MPFR_RNDZ,
    404                         MPC_PREC_RE (rop_sin)
    405                         + (MPC_RND_RE (rnd_sin) == MPFR_RNDN));
    406 
    407             if (ok) {
    408                /* imaginary part of sine */
    409                mpfr_mul (csh, c, sh, MPFR_RNDN);
    410                ok = (!mpfr_number_p (csh))
    411                      || mpfr_can_round (csh, prec - 2, MPFR_RNDN, MPFR_RNDZ,
    412                            MPC_PREC_IM (rop_sin)
    413                            + (MPC_RND_IM (rnd_sin) == MPFR_RNDN));
    414             }
    415          }
    416 
    417          if (rop_cos != NULL && ok) {
    418             /* real part of cosine */
    419             mpfr_mul (c, c, ch, MPFR_RNDN);
    420             ok = (!mpfr_number_p (c))
    421                   || mpfr_can_round (c, prec - 2, MPFR_RNDN, MPFR_RNDZ,
    422                         MPC_PREC_RE (rop_cos)
    423                         + (MPC_RND_RE (rnd_cos) == MPFR_RNDN));
    424 
    425             if (ok) {
    426                /* imaginary part of cosine */
    427                mpfr_mul (s, s, sh, MPFR_RNDN);
    428                mpfr_neg (s, s, MPFR_RNDN);
    429                ok = (!mpfr_number_p (s))
    430                      || mpfr_can_round (s, prec - 2, MPFR_RNDN, MPFR_RNDZ,
    431                            MPC_PREC_IM (rop_cos)
    432                            + (MPC_RND_IM (rnd_cos) == MPFR_RNDN));
    433             }
    434          }
    435       } while (ok == 0);
    436 
    437       if (rop_sin != NULL) {
    438          inex_re = mpfr_set (mpc_realref (rop_sin), sch, MPC_RND_RE (rnd_sin));
    439          if (mpfr_inf_p (sch))
    440            inex_re = mpc_fix_inf (mpc_realref (rop_sin), MPC_RND_RE (rnd_sin));
    441          inex_im = mpfr_set (mpc_imagref (rop_sin), csh, MPC_RND_IM (rnd_sin));
    442          if (mpfr_inf_p (csh))
    443            inex_im = mpc_fix_inf (mpc_imagref (rop_sin), MPC_RND_IM (rnd_sin));
    444          inex_sin = MPC_INEX (inex_re, inex_im);
    445       }
    446       else
    447          inex_sin = MPC_INEX (0,0); /* return exact if not computed */
    448 
    449       if (rop_cos != NULL) {
    450          inex_re = mpfr_set (mpc_realref (rop_cos), c, MPC_RND_RE (rnd_cos));
    451          if (mpfr_inf_p (c))
    452            inex_re = mpc_fix_inf (mpc_realref (rop_cos), MPC_RND_RE (rnd_cos));
    453          inex_im = mpfr_set (mpc_imagref (rop_cos), s, MPC_RND_IM (rnd_cos));
    454          if (mpfr_inf_p (s))
    455            inex_im = mpc_fix_inf (mpc_imagref (rop_cos), MPC_RND_IM (rnd_cos));
    456          inex_cos = MPC_INEX (inex_re, inex_im);
    457       }
    458       else
    459          inex_cos = MPC_INEX (0,0); /* return exact if not computed */
    460 
    461       mpfr_clear (s);
    462       mpfr_clear (c);
    463       mpfr_clear (sh);
    464       mpfr_clear (ch);
    465       mpfr_clear (sch);
    466       mpfr_clear (csh);
    467 
    468       /* restore the exponent range, and check the range of results */
    469       mpfr_set_emin (saved_emin);
    470       mpfr_set_emax (saved_emax);
    471       if (rop_sin != NULL)
    472         {
    473           inex_re = mpfr_check_range (mpc_realref (rop_sin),
    474                                       MPC_INEX_RE (inex_sin),
    475                                       MPC_RND_RE (rnd_sin));
    476           inex_im = mpfr_check_range (mpc_imagref (rop_sin),
    477                                       MPC_INEX_IM (inex_sin),
    478                                       MPC_RND_IM (rnd_sin));
    479           inex_sin = MPC_INEX (inex_re, inex_im);
    480         }
    481       if (rop_cos != NULL)
    482         {
    483           inex_re = mpfr_check_range (mpc_realref (rop_cos),
    484                                       MPC_INEX_RE (inex_cos),
    485                                       MPC_RND_RE (rnd_cos));
    486           inex_im = mpfr_check_range (mpc_imagref (rop_cos),
    487                                       MPC_INEX_IM (inex_cos),
    488                                       MPC_RND_IM (rnd_cos));
    489           inex_cos = MPC_INEX (inex_re, inex_im);
    490         }
    491 
    492       return (MPC_INEX12 (inex_sin, inex_cos));
    493    }
    494 }
    495